

#ifndef _CDS_LD_H_
#define _CDS_LD_H_

#include "Tuna.hpp"
#include "Equations/ScalarEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class CDS_LD;

  template <typename T>
  struct Typeinfo <CDS_LD<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <CDS_LD<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <CDS_LD<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };
  
/*!
 ***************************************************************************
 * Adaptor for central differences for energy equation.
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Mon Mar 17 19:55:58 GMT 2008 ]
 ***************************************************************************
 */
template<typename Tprec, int Dim>
class CDS_LD : public ScalarEquation< CDS_LD< Tprec, Dim > >
{
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aE;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aW;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aN;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aS;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aF;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aB;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::aP;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::sp;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::dx;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::dy;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::dz;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::bi;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::ei;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::bj;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::ej;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::bk;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::ek;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::dt;    
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::applyBoundaryConditions1D;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::applyBoundaryConditions2D;
    using GeneralEquation<ScalarEquation<CDS_LD<Tprec, Dim > > >::applyBoundaryConditions3D;
    
    using ScalarEquation<CDS_LD<Tprec, Dim > >::Gamma;
    using ScalarEquation<CDS_LD<Tprec, Dim > >::phi_0;
    using ScalarEquation<CDS_LD<Tprec, Dim > >::u;
    using ScalarEquation<CDS_LD<Tprec, Dim > >::v;
    using ScalarEquation<CDS_LD<Tprec, Dim > >::w;
    
public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim>::huge ScalarField;

    CDS_LD() : ScalarEquation<CDS_LD<prec_t, Dim> >() { }
    ~CDS_LD() { }; 

    inline bool calcCoefficients1D() { };
    inline bool calcCoefficients2D();
    inline bool calcCoefficients3D() { };
    inline void printInfo() { std::cout << " CDS_LD "; }
};


template<typename Tprec, int Dim>
inline bool CDS_LD<Tprec, Dim>::calcCoefficients2D () 
{
    prec_t Gdy_dx = Gamma * dy / dx;
    prec_t Gdx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i) 
	for (int j = bj; j <= ej; ++j) {
	    aE (i,j) = Gdy_dx; 
    	    aW (i,j) = Gdy_dx; 
	    aN (i,j) = Gdx_dy; 
	    aS (i,j) = Gdx_dy; 
	    aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt;
	    sp (i,j) = phi_0(i,j) * dxy_dt - u(i,j) * dx * dy;    
	}
    applyBoundaryConditions2D();
    return 0;
}

} // Tuna namespace


#endif //_CDS_LD_H_

















